The ‘anglr’ package illustrates some generalizations of GIS-y and topology tasks in R with “tables”. See the package project for more information.
Get some maps and plot in 3D - in plane view, or globe view.
NOTE: these plots are interactive, but performance and quality will be better (for now) if run locally.
Set up a simple world countries data set.
library(rgl)
simpleworld <- rnaturalearth::countries110
library(raster)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
## convert to triangles and plot
library(anglr)
#>
#> Attaching package: 'anglr'
#> The following object is masked _by_ '.GlobalEnv':
#>
#> simpleworldYou must enable Javascript to view this page properly.
rgl::rgl.clear()
sids <- sf::read_sf(system.file("shape/nc.shp", package="sf"))
ex <- extent(sf::st_bbox(sids)[c(1, 3, 2, 4)]) + 5
gl <- graticule::graticule(seq(xmin(ex), xmax(ex), length = 15),
seq(ymin(ex), ymax(ex), length = 8))
## convert to triangles, but wrap onto globe then plot
smesh <- anglr(sids)
plot(globe(smesh))
mgl <- anglr(gl)
mgl$o$color_ <- "black"
plot(globe(mgl), lwd = 2)
rgl::rglwidget()You must enable Javascript to view this page properly.
It’s trivial to have holes, because there are no holes, we have a true surface, composed of 2D primitives (triangles).
rgl::rgl.clear()
library(spbabel)
data(holey)
## SpatialPolygonsDataFrame
sph <- sp(holey)
glh <- anglr(sph)
plot(glh)
rgl::rglwidget()You must enable Javascript to view this page properly.
rgl::rgl.clear()
linehouse <- as(sph, "SpatialLinesDataFrame")
plot(anglr(linehouse))
#> Joining, by = "object_"
#> Joining, by = "segment_"
rgl::rglwidget()You must enable Javascript to view this page properly.
rgl::rgl.clear()
lmesh <- anglr(as(simpleworld, "SpatialLinesDataFrame"))
plot(globe(lmesh))
#> Joining, by = "object_"
#> Joining, by = "segment_"
rgl::rglwidget()You must enable Javascript to view this page properly.
Rgl mesh3d objects that use “triangle” primitives are supported.
rgl::rgl.clear()
dod <- anglr(dodecahedron3d(col = "cyan"))
octo <- anglr(translate3d(octahedron3d(col = "blue"), 6, 0, 0))
plot(dod, col = viridis::viridis(5)[1], alpha = 0.3)
plot(octo, col = viridis::viridis(5)[5], alpha = 0.3)
bg3d("grey")
rgl::rglwidget()You must enable Javascript to view this page properly.
To complete the support for these rgl objects we need quads, and to allows a mix of quads and triangles in one data set (that’s what extrude3d uses). Extrusion is a bit limiting so it’s unclear how useful this is (yes it is common, though). See rgl::extrude3d for the most readily available method in R.
And points work! (Don’t laugh).
rgl::rgl.clear()
library(anglr)
mpts <- as(as(simpleworld, "SpatialLinesDataFrame"), "SpatialMultiPointsDataFrame")
plot(anglr(mpts))
#> Joining, by = "object_"
#> Joining, by = "branch_"
rgl::view3d(theta = 25, phi = 3)
rgl::rglwidget()You must enable Javascript to view this page properly.
The trip package includes a ‘walrus818’ data set courtesy of Anthony Fischbach. Zoom around and see if you can find them.
Here we also use the finite element aspect of our topology representation. We can forget about complicated algorithms to cut and splice lines and just throw away triangles south of a chosen latitude. (Still needs work, but we’ve proved it works well enough. )
rgl::rgl.clear()
library(trip)
library(anglr)
data(walrus818)
library(graticule)
prj <-"+proj=laea +lon_0=0 +lat_0=90 +ellps=WGS84"
gr <- graticule(lats = seq(40, 85, by = 5), ylim = c(35, 89.5), proj = prj)
w <- spTransform(subset(simpleworld, coordinates(simpleworld)[,2] > 10), prj)
walrus <- spTransform(walrus818, prj)
gr$color_ <- "black"
rgl::par3d(windowRect = c(100, 100, 912 + 100, 912 +100))
plot(anglr(gr))
#> Joining, by = "object_"
#> Joining, by = "segment_"
w$color_ <- sample(viridis::inferno(nrow(w)))
wmap <- anglr(w)
tXv <- wmap$tXv %>% dplyr::inner_join(wmap$v)
#> Joining, by = "vertex_"
keep <- tibble::tibble(triangle_ = unique(tXv$triangle_[rgdal::project(as.matrix(tXv[c("x_", "y_")]), prj, inv = TRUE)[,2] > 35]))
wmap$t <- wmap$t %>% dplyr::inner_join(keep)
#> Joining, by = "triangle_"
## trim by primitives
## we mashed the graph and bit but it works, later we
## can propagate the filter ...
plot(wmap, specular = "black")
plot(anglr(walrus))
#> Joining, by = "object_"
#> Joining, by = "segment_"
um <- structure(c(0.934230506420135, 0.343760699033737, 0.0950899347662926,
0, -0.302941381931305, 0.905495941638947, -0.297159105539322,
0, -0.188255190849304, 0.24880850315094, 0.950081348419189, 0,
0, 0, 0, 1), .Dim = c(4L, 4L))
par3d(userMatrix = um)
rgl::rglwidget()You must enable Javascript to view this page properly.
A polygon layer may be overlaid with a raster layer using the z argument.
rgl::rgl.clear()
topo <- raster::raster(system.file("extdata", "gebco1.tif", package = "anglr"))
poly <- subset(simpleworld, sovereignt == "India")
prj <- "+proj=laea +lon_0=80 +lat_0=45 +datum=WGS84"
poly_z <- anglr(poly, z = topo * 68, max_area = 0.02)
#> Warning in anglr.PATH(PATH(x), z = z, ..., type = type, max_area =
#> max_area):
poly_z$v[c("x_", "y_")] <- rgdal::project(as.matrix(poly_z$v[c("x_", "y_")]), prj)
poly_z$meta$proj <- prj
plot(poly_z)
rgl::rglwidget()